library(tidyverse)
install.packages("cowplot")
library(cowplot)
library(ggpubr)
library(pheatmap)
library(RColorBrewer)
library(vegan)
install.packages("reshape2")
library(reshape2)
install.packages("corrplot")
library(corrplot)
library(kableExtra)
library(agricolae)
install.packages("tidyr")
library(tidyr)
library(magrittr)
library(ggplot2)
library(readr)
library(dplyr)
#install.packages("tibble")
library(tibble)
# 设置默认主题
theme_set(theme_bw())

#biolog <- read_csv("ratio/data/biolog.csv")
biolog <- read_csv("C:/Users/兰兰子/Downloads/ratio-master/data/biolog.csv")

#qPCR_data <- read_csv("ratio/data/qPCR.csv")
qPCR <- read_csv("C:/Users/兰兰子/Downloads/ratio-master/data/qPCR.csv")

#mono_data <- read_csv("ratio/data/mono.csv")
mono <- read_csv("C:/Users/兰兰子/Downloads/ratio-master/data/mono.csv")
#head(biolog)
head(biolog)
head(qPCR)
head(mono)
carbon <- read_csv("C:/Users/兰兰子/Downloads/ratio-master/data/carbon.csv")
head(carbon)

#######################12.5################

mono_data <- mono %>% 
  melt(
    id.vars=c("plate","carbon_id"),
    variable.name ="Target.Name",
    value.name="Quantity_mono"
    )

cocu_data <- qPCR %>% 
  select(
    plate,carbon_id,ratio0,EC,PP
    ) %>% 
  melt(
    id.vars=c("plate","carbon_id","ratio0"),
    variable.name ="Target.Name",
    value.name="Quantity_cocu"
    )

data_all <- merge(mono_data, cocu_data, by = c("carbon_id","Target.Name","plate"),all=T) %>% 
  filter(carbon_id!="1")



qPCR_data <- qPCR %>%
  mutate(ratio0 = factor(ratio0, levels = c("less","equal","more")))

# 标准化
biolog_24h <- biolog %>% 
  mutate(ratio0 = factor(ratio0, levels = c("none","less","equal","more","all"))) %>%
  group_by(plate,ratio0) %>% 
  mutate(A590=A590-A590[carbon_id==1],A750=A750-A750[carbon_id==1]) %>%   # 将阴性对照设为零
  filter(carbon_id!=1) %>%
  ungroup()

biolog_mono_24h <- biolog_24h %>% 
  filter(ratio0 %in% c("none","all")) %>% 
  mutate(species=factor(ratio0,levels = c("all","none"),labels = c("E. coli","P. putida"))) %>% 
  dplyr::select(-ratio0)

biolog_coculture_24h <- biolog_24h %>% 
  filter(ratio0 %in% c("less","equal","more")) %>%
  mutate(ratio0 = factor(ratio0, levels = c("less","equal","more")))


#12.6.1
M_A590_24h <- biolog_24h %>% mutate(sample=paste(ratio0,plate,sep="-")) %>%
  dplyr::select(sample,carbon_id,A590) %>%
  spread(key=sample,value=A590) %>%
  as.data.frame() %>%
  tibble::column_to_rownames(var="carbon_id")

k3 <- cutree(hclust(dist(M_A590_24h)),k=3)

carbon_group <-  data.frame(usage=k3) %>%
  rownames_to_column(var="carbon_id") %>%
  mutate(carbon_id=as.numeric(carbon_id)) %>%
  mutate(usage=paste("U",usage,sep=""))

carbon_name <- left_join(carbon, carbon_group)

#12.6.2

biolog_mono_A590_24h <- biolog_mono_24h %>% 
  dplyr::select(plate,carbon_id,species,A590) %>% 
  spread(species,A590) 

PP_prefered <- biolog_mono_A590_24h %>% 
  group_by(carbon_id) %>%  
  summarise(p=t.test(`P. putida`,`E. coli`,alternative = "greater")$p.value) %>% 
  filter(p<0.05)

EC_prefered <- biolog_mono_A590_24h %>% 
  group_by(carbon_id) %>%  
  summarise(p=t.test(`P. putida`,`E. coli`,alternative = "less")$p.value) %>% 
  filter(p<0.05)

carbon_prefer <- data.frame("carbon_id"=carbon_name$carbon_id,
                            "prefer"="None",
                            stringsAsFactors = F)

carbon_prefer[carbon_prefer$carbon_id %in% EC_prefered$carbon_id,"prefer"] <- "EC"

carbon_prefer[carbon_prefer$carbon_id %in% PP_prefered$carbon_id,"prefer"] <- "PP"

carbon_name <- left_join(carbon_name, carbon_prefer)

#12.6.3
carbon_name %>% 
  left_join(carbon_prefer) |> 
  kableExtra::kable()

#12.7(1图)
ratio1 <- qPCR_data %>% filter(ratio0 %in% c("less","equal","more")) %>%
  complete(ratio0,carbon_id,plate) %>% 
  group_by(ratio0,carbon_id) %>% 
  dplyr::select(ratio0,plate,carbon_id,ratio1) %>% 
  mutate(ratio1_mean=mean(ratio1,na.rm = T)) %>% 
  mutate(ratio1=ifelse(is.na(ratio1),ratio1_mean,ratio1)) %>% 
  dplyr::select(-ratio1_mean)

mono_A590 <- biolog_mono_24h %>% 
  group_by(carbon_id,species) %>% 
  summarise(A590=mean(A590)) %>% 
  spread(key="species",value="A590") 

A590_caculated <- left_join(ratio1,mono_A590) %>% 
  mutate(A590_cac=(`P. putida`+ratio1*`E. coli`)/(1+ratio1))

social <- biolog_coculture_24h %>% 
  dplyr::select(plate,carbon_id,ratio0,A590) %>%
  left_join(A590_caculated) %>%
  group_by(carbon_id,ratio0) %>% 
  mutate(p_pos=t.test(x=A590,y=A590_cac,alternative = "greater")$p.value,
         p_neg=t.test(x=A590,y=A590_cac,alternative = "less")$p.value) %>%
  mutate(social_type=ifelse(
    p_pos<0.05,"+",
    ifelse(p_neg<0.05,"-","unresolved"))
  ) %>% 
  ungroup() %>%
  dplyr::select(carbon_id,ratio0,social_type,p_pos,p_neg) %>%
  unique() %>%
  mutate(ratio0=factor(ratio0,levels = c("less","equal","more")))

table(social$social_type) %>% barplot(col=c("blue","red","grey"))
#12.7(2图)

social_qpcr <- data_all %>% 
  group_by(carbon_id,Target.Name,ratio0) %>% 
  summarise(p_pos=t.test(x=log10(Quantity_cocu),y=log10(Quantity_mono),alternative = "greater")$p.value,
            p_neg=t.test(x=log10(Quantity_cocu),y=log10(Quantity_mono),alternative = "less")$p.value) %>%
  mutate(social_type=ifelse(
    p_pos<0.05,"+",
    ifelse(p_neg<0.05,"-","unresolved"))
  ) %>% ungroup() %>%
  mutate(ratio0=factor(ratio0,levels = c("less","equal","more")))

table(social_qpcr$social_type) %>% barplot(col=c("blue","red","grey"))

#12.8
merged <- left_join(biolog_coculture_24h,qPCR_data) %>% 
  left_join(social) %>% 
  left_join(carbon_name) %>%
  filter(!is.na(ratio1))

social_qpcr <- social_qpcr %>% select(social_type,carbon_id,ratio0)

merged_qpcr <- qPCR_data %>% 
  left_join(social_qpcr) %>% 
  left_join(carbon_name) %>%
  filter(!is.na(ratio1)) %>% 
  mutate(social_type=factor(social_type,levels = c('unresolved','+','-'))) %>% 
  mutate(prefer=factor(prefer,levels = c('none','EC','PP')))

merged_qpcr$EC <- log10(merged_qpcr$EC)
merged_qpcr$PP <- log10(merged_qpcr$PP)

#12.8.1
merged <- merged %>% filter(ratio0 %in% c("less","equal","more"))

par(mfrow=c(3,4))

hist(merged$EC)
qqnorm(merged$EC)
hist(log10(merged$EC))
qqnorm(log10(merged$EC))

hist(merged$PP)
qqnorm(merged$PP)
hist(log10(merged$PP))
qqnorm(log10(merged$PP))

hist(merged$ratio1)
qqnorm(merged$ratio1)
hist(log10(merged$ratio1))
qqnorm(log10(merged$ratio1))

